Eigenstate Estimation for the Bardeen-Cooper-Schrieffer (BCS) Hamiltonian 
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We show how multi-level BCS Hamiltonians of finite systems in the strong pairing interaction 
regime can be accurately approximated using multi-dimensional shifted harmonic oscillator Hamil- 
tonians. In the Shifted Harmonic Approximation (SHA), discrete quantum state variables are ap- 
proximated as continuous ones and algebraic Hamiltonians are replaced by differential operators. 
Using the SHA, the results of the BCS theory, such as the gap equations, can be easily derived with- 
out the BCS approximation. In addition, the SHA preserves the symmetries associated with the 
BCS Hamiltonians. Lastly, for all interaction strengths, the SHA can be used to identify the most 
important basis states - allowing accurate computation of low-lying eigenstates by diagonalizing 
BCS Hamiltonians in small subspaces of what may otherwise be vastly larger Hilbert spaces. 
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The traditional method of finding eigenvalues of a 
Hamiltonian H({X^}) (expressed as a polynomial in the 
elements {X„} of a Lie algebra g) is by diagonalization. 
However, in realistic many-body systems the Hamilto- 
nian matrices can be huge. The problem is then to find 
an approximation such that the salient features of the 
Hamiltonian are retained. In this letter, the so-called 
Shifted Harmonic Approximation (SHA), introduced by 
Chen et al. is developed and extended to many de- 
grees of freedom. The key principle behind the SHA is 
to replace discrete quantum state variables by contin- 
uous ones. Algebraic Hamiltonians are then replaced 
by differential operators. This approach offers new in- 
sights even for well studied systems such as those with 
a Bardccn-Coopcr-Schricffcr (BCS) Hamiltonian 0, 
which in general cannot be solved exactly. The tradi- 
tional BCS approximation provides accurate results in 
the thermodynamic limit but violates particle-number 
conservation. For finite systems, this is a major source 
of inaccuracies but its effects can be reduced by number 
conserving extensions of the BCS theory such as 0, [1[ ■ 

Recent studies of superconductivity in metallic nano- 
grams [(| and atomic nuclei Q have led to a revival of in- 
terest in the Richardson- Gaudin approach [1, @. Classes 
of BCS Hamiltonians with level-independent interactions 
are shown to be integrable and solvable by means of an 
algebraic Bethe ansatz. However, the numerical solutions 
are challenging to compute jic| and the eigenstates are 
not easy to use. Moreover, among the set of BCS Hamil- 
tonians, there is only a small number of special cases fill ] 
that are solvable by the Richardson-Gaudin method. 

Here, using the SHA which is number conserving, we 
show that a general fc-level BCS Hamiltonian can be ap- 
proximated as a (k — l)-dimcnsional shifted oscillator 



Hamiltonian. Accurate approximations of the low-lying 
eigenstates are then easily obtained in the strong interac- 
tion regime. In the weak interaction regime, the SHA can 
also be used to identify the most important basis states 
for computing the low-lying eigenstates accurately. 

Consider an irreducible representation (irrep) of the 
su(2) algebra on the Hilbert space spanned by basis states 
{|m),m = —j, Any state \4>) in this Hilbert 

space, e.g., an eigenstate of a Hamiltonian in the su(2) 
algebra, can be expressed as a linear combination of the 
basis states \<f>) = J2 m l m )( m l^) = Em l m )<M m )j where 
the coefficient <f>(m) = (m\4>) is a discrete distribution of 
m. The action of the su(2) operators on such a distribu- 
tion, defined by Jk4>{m) = (m\Jk\<j>), is then 

J z (j){m) = m0(m), (1) 
J±4>{m) = ^{jTm + l){j±m)4>{m T 1). (2) 

For large values of j and for a state for which </>(m) varies 
slowly with the discrete variable m, we can now make 
the continuous variable approximation of extending m to 
continuous values and replacing <p(m) by a smooth func- 
tion ijj(x), defined such that ip( x ) = 4>{m) when x — m/j. 
We can then use the identity ip(x^f -.) = exp ( = Fy^)V'( a; ) 
and, assuming the expansion of exp (-fj^:)'4'{ x ) to be 
rapidly convergent, make the approximation 
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Note that we have omitted the 1/j term to obtain 
— x 2 in cq. ([3]) . This term is negligible for large values 
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of j, but could be included in a more complete calcula- 
tion. If the function ijj(x) is (i) slow varying, (ii) local- 
ized about a value x D and (iii) vanishes when |a:| — > 1, we 
can make the Shifted Harmonic Approximation (SHA). 
In this approximation, the action of an su(2) operator, 
such as J± in eq. (O, on i\){x) is obtained by expanding it 
about x up to bilinear terms. Similarly, a Hamiltonian 
that is quadratic in the elements of an su(2) algebra and 
has low-lying eigenfunctions that satisfy the SHA crite- 
ria can be mapped to a harmonic oscillator Hamiltonian 
'Hsha that is bilinear in (x — x ) and d/dx. 

Now consider a multi-level BCS Hamiltonian consist- 
ing of fermions in k single-particle energy levels. The op- 
erators (fl/jp) create (annihilate) a fermion in a state 
/i p at level p, and the operators for the corresponding 
time-reversed states are ajj p (a^ p ). As shown by Kerman 
et. al. 0] , these operators can be combined to form su(2) 
quasi-spin operators 

^ = IE«s-^U Jl = \Y.«- ( 4 ) 

Mp>0 Hp 

Here, the operator creates a pair of particles in time- 
reversed states at level p. Together with the pair annihi- 
lation operators, = (</+), these quasi-spin operators 
belong to an su(2) su(2) ® . . . algebra with commuta- 
tion relations J q _] = 2JPS pq and [J£,J£] = ±J^S pq . 
In this formalism, the BCS Hamiltonian is written as 



P =i 



k 

p,q 



(5) 



where n p = 2(Jf + j p ) is the particle number operator 
for the level with single particle energy e p . The operator 
J+J'L scatters a pair of particles from level q to level p 
and G pq is the corresponding interaction strength. 

The Hamiltonian (0 conserves both particle number 
and the number of paired particles. Without loss of gen- 
erality, we consider a system with no unpaired particles. 
For level p, let {\j P m p }} denote the basis states for the 
irreducible su(2) p representation for which \j p , m p = —j p ) 
is the zero-pair state and m p increases by one with every 
added pair to reach the value m p = j p , when the level 
is completely filled. Basis states for the fc-level pairing 
model with no unpaired particles, are then defined by 

l m ) = \ji m i) ® \j2m2) <& • • • ® lifemfe). 

Eigcnstatcs |<I> 1 ) of the pairing Hamiltonian ([5]), with a 
fixed pair number A, are given by linear combinations of 
the basis states, |m), for which 5Z p =i 

m p = N- \ A max 

where A max is the maximum number of pairs possible in 
the system. If we plot the set of allowed basis states |m) 
for the A-pair system as points on an {mi, 7712, . . . , mk} 
grid, these points lie on a (k— l)-dimcnsional hypcrplane. 
The direction orthogonal to this hypcrplane is described 
as spurious because there is no dynamics associated with 
it when A is fixed. See FIG.Q]for a sample 2- level system. 




FIG. 1: (Color online) A sample two-level system. The com- 
ponents of the exact ground state eigenvectors ($ 1 (mi, 7712)) 
are indicated with dots on an (mi, 7712) plane for (a) large 
and (b) small interactions. The solid and dash lines are the 
corresponding SHA eigenfunctions. The ground state for a 
different N is shown by (c). The directions of the transformed 
coordinates, £1 and £ sp = £2, with their origins at the point 
O ~ (jix \, J2X02) are indicated for (a). 



To apply the SHA to the pairing Hamiltonian we 
define 'H$ i (m) = (m\H\& 1 ). Then, assuming that the 
low-lying eigenfunctions $' (m) vary slowly with m, we 
define k continuous variables x p — m p / j p , and the contin- 
uous eigenfunction ^'(x) as before. Assuming also that 
the wave functions, ^(x), are localized around a point 
x D , and that the criteria for the validity of the SHA, listed 
above, are satisfied, we expand the Hamiltonian % up to 
bilinear terms in x' p = (x p — x op ) and -^p- to obtain 
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which is essentially the Hamiltonian for a coupled k- 
dimensional harmonic oscillator with the origin shifted 
to x op . It contains: an inverse mass tensor A pq , a spring 
constant tensor B pq , a set of shifts D p and a constant E. 
Their components, defined in terms of K r = (1 — x 2 or ) and 
T pq = G pq j p j qy /K p ~K q ~, are given by 
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(8) 

(9) 
(10) 
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The conservation of particle number in the SHA formal- 
ism is verified by showing that the SHA representation of 



3 



the number operator n = 2^2 p {J% + j p ) commutes with 
Usha- Thus, it is appropriate to make a change of vari- 
ables, jpx' p — > £i, such that (denoted £ sp in FIG. [1]) 
is an iV-dependent constant and the Hamiltonian ([5]) be- 
comes that of a system of (k — 1) harmonic oscillators, 

Observe that while the dynamics of the oscillator 
is (k — l)-dimensional, the tensors A and B are k- 
dimensional. However, the inverse mass tensor A is de- 
termined to have an eigenvector parallel to the spurious 
direction with zero eigenvalue; this implies that the cor- 
responding vibrational mass is infinite, consistent with 
N being a good quantum number in the SHA. To eval- 
uate the quantities in eqs. ([71 - 110[) . we need the numer- 
ical values of x op . The point x Q is naturally defined as 
the minimum of the harmonic oscillator potential on the 
(k — l)-dimensional hyperplane (see FIG. [1}. Once the 
pair number N is selected, the numerical values of x op 
are determined by the D p shift functions of Eq. ((9|) . 

Having determined x op , we can use the transformation 
that diagonalizcs A to obtain the dynamics on the hyper- 
plane. The properties of the (k— l)-dimcnsional oscillator 
on the hyperplane are solved using normal-mode theory. 
The eigenvalues of the SHA Hamiltonian are 

fc-i 

S v = ^2(y i + ^)u i + E, (12) 
j=i 

where w, = and v = {vi} is a set of integers 

indicating the number of oscillator quanta in mode i. 
The corresponding set of SHA eigenfunctions are 

fc-i 

MO = vU {r^ lVl \^y h H v ^yr^ ? (13) 
t=i 

where H Vi are Hermite polynomials, (7; = are the 

SHA widths, and rj is a normalization factor. From here, 
we can approximate the coefficients $ J (m) = (m|$ J ), 
and hence the eigenfunction, in the original discrete ba- 
sis by evaluating ^ w at the points £ corresponding to 
m. We refer to these approximate eigenfunctions of the 
Hamiltonian H, as the SHA basis. 

It is worth noting that the quantity \{x op + 1) in the 
SHA can be interpreted as the mean fractional occupancy 
of a single-particle level p in parallel with v p in BCS the- 
ory Q. Similarly, the SHA shift equations for x op cor- 
respond to the BCS gap equations for v p . In addition, 
the SHA energy E is almost identical to the BCS ground 
state energy. It has been shown, in several model calcula- 
tions, that including higher order corrections in i lowers 
the SHA ground-state energy in the strong interaction 
regime below that of the BCS approximation. Thus, we 



obtain an insightful interpretation of the SHA treatment 
of the pairing model as an extension of the BCS method 
to a number conserving approximation which takes ac- 
count of the fluctuations of the particle number in each 
single-particle level about its mean BCS value. 

The transformed ^-coordinates for a 2-level model are 
illustrated in FIG.[TJ The SHA eigenfunction (line) corre- 
sponding to a large interaction (compared to single parti- 
cle energies spacing) , indicated as (a) , are in good agree- 
ment with the exact components of the eigenvectors given 
by diagonalization. Similar accuracy is obtained for the 
next few higher-energy states (not shown) . If greater pre- 
cision is required, an even more accurate description of 
the low-lying eigenstates can be obtained by diagonal- 
izing the BCS Hamiltonian in a subspace spanned by a 
small number of SHA basis states. For weaker interac- 
tions, as in (b), the SHA does not predict the components 
of the eigenvectors accurately as in (a). This is the regime 
in which the conditions for the validity of the SHA are 
not well satisfied. Nevertheless, some SHA predictions, 
such and <Ji, remain accurate - a subtlety not yet 

fully understood. Thus, we can use these predictions 
to identify a small subset of basis states that contribute 
significantly to the low-lying eigenstates in the weak in- 
teraction regime and also obtain very accurate results for 
them by diagonalizing small Hamiltonian matrices. 

To illustrate the effectiveness of the SHA, we consider 
a system with four degenerate single-particle energy lev- 
els. This relatively small system is selected so that exact 
eigenstates can be obtained by diagonalization. Based on 
other applications of the SHA [12j, we expect the SHA 
to be even more accurate and effective in application to 
systems of single-particle levels of higher multiplicities. 

Exact and SHA-estimated excitation energies of a sam- 
ple 4-lcvel model with N = 28 pairs, j = [7, 8, 9, 10] and 
e = [0.5,2.3,6.1,7.3] arc shown in FIG. [2 A simple ar- 
bitrary rule for G is used: G pq = (2.0 — 0.1 \e p — e q \)g, 
where g controls the interaction strength. Exact results 
are obtained by diagonalizing the 3231 x 3231 Hamil- 
tonian matrix. For this 4-level model, the SHA oscil- 
lator in the strong interaction regime is 3-dimensional. 
The number of oscillator quanta in each mode is given 
by {ni,n 2 ,n 3 }. The excitation energy A£{„ lin2jn3 } = 
£{ni.n 2 ,ri3} — £{o,o,o} f° r the low-lying states are shown. 
From the figure, we see that the SHA-predicted excita- 
tion energies ('+', 'x') are in good agreement with the 
exact results from diagonalization ('©') for a wide range 
of interactions. Note that while the SHA excitation en- 
ergies for g ~ 0.05 are less accurate than for other values 
of g, the trends in how the excitation energies vary with 
interaction strength are still closely captured by the SHA. 

Lastly, we show in TABLE U the low-lying excitation 
energies obtained from the SHA and by diagonalizing the 
BCS Hamiltonian in the space spanned by the most im- 
portant basis states identified by the SHA. For points (1) 
and (2) in FIG.HJ the lowest 50 and 300 su(2) basis states 
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FIG. 2: (Color online) The SHA predicted excitation energies 
for {m, 7i2, ri3}={l, 0, 0}, {0, 1, 0} and {0,0,1} are indicated 
with ' + ' and those for {2,0,0} are indicated with 'x'. The 
exact lowest four excitation energies are given by 'o'. We 
mark three points of interest (1) g = 0.010, (2) g = 0.045 and 
(3) g — 0.150. Note that the SHA excitation energies corre- 
sponding to {2, 0, 0} cross-over those of {0, 0, 1} at g ~ 0.05, 
indicating the system transiting from one regime to another. 



TABLE I: The lowest few excitation energies for (1) g = 0.010 
(2) g = 0.045 and (2) g = 0.150 as predicted by the SHA and 
diagonalizing in the most relevant subspaces (indicated by 
Diag) are shown. The 'Diag' results are given up to the digit 
that agrees with the corresponding exact result. All SHA 
predictions are given to two decimal places in ( ) with the 
number of oscillator quanta given in { } in the subscript. 



(1) 3 = 0.010 
Diag S u(2) (SHA) 



(2) g = 0.045 
Diag.u(2) (SHA) 



(3) g = 0.150 
L>ia flS HA(SHA) 



3.650313 (3.61) { i i0 ,o} 3.60 (3.07) { i i0 ,o> 15.03 (15.06) { i i0 ,o} 

6.9484 (6.90) {0 ,i,o } 5.258 (5.19) {0 ,i, } 16.18 (16.22) {0 ,i,o } 

7.38050 (7.21) {2 ,o, 0} 7.65 (6.14) {2 , , 0} 18.1 (18.08) {0 ,o,i } 

9.4210 (9.38) {0 ,o,i> 7.8 (7.20) {0 , , 1} 29.5 (30.12) {2 , , 0} 

10.5531 (10.51){i,i, 0} 8.37 (8.26) {M , 0} 30.7 (31.28) {M , 0} 



are used respectively. For point (3), we used the lowest 
286 SHA wave functions as a basis. 

The results obtained, cf. last column of TABLE[H show 
the SHA to be very successful for deriving low-energy 
spectra of BCS Hamiltonians in the strong interaction 
regime in which it is most valid. They also show the 
SHA to be a good first-order approximation in general. 
As TABLE [T] indicates, accurate results can be obtained 
for any interaction by using the SHA to select relatively 
small subsets of basis states for diagonalizations. 

The SHA predicts essentially the same mean level oc- 



cupancies in the ground state as the BCS approximation. 
In addition, it gives the fluctuations in these occupan- 
cies in a manner that conserves particle number. It also 
conserves the su(2) <E> su(2) . . . symmetry of the pairing 
Hamiltonian (defined by the values of the j p quantum 
numbers). The SHA gives the low-energy states of all 
irreps of this symmetry group. This is in contrast to 
the BCS approximation which is only designed to give 
an approximation for the ground state and quasi-particle 
approximations for the low-energy states of neighbouring 
odd-particle systems. States of maximal su(2) ® su(2) . . . 
symmetry are unbroken-pair states, whereas the broken- 
pair states of other irreps have unpaired particles in one 
or more single-particle levels. This reduces the number 
of states available to the paired particles so that these ir- 
reps are obtained by reducing the quasi-spin of each level 
p by the replacement j p — > j p — \ for each unpaired par- 
ticle in the level. The states of such irreps are handled in 
the same way in the SHA, except for the different values 
of the quasi-spins. To conclude, we note the significant 
possibility that the continuous variable approximation, 
underlying the SHA, has the potential to be applied to 
derive other solvable differential equations. This poten- 
tial remains to be explored. 
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